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(54) Titie: MULTIEXPONENTIAL SIGNAL PROCESSING METHOD AND APPARATUS 
(57) Abstract 

For signal processing of tran- 
sients such as multiexponential decays, 
a transf rm operator (e.g., matrix) is 
constructed by emphasizing linear res- 
lution. rather than using fitting rou- 
tines that attempt to find an estimate of 
an imknown model that fits data. Use 
of the transform operator to process 
multiexponential signals produces out- 
puts that are a better estimate of the 
unknown model or of some segment 
of the unknown model. 
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Title: 



Multiexponential Signal Processing 
Method and Apparatus 



SPECIFICATION 



CROSS-REFERENCE TO RELATED APPLICATION 

This application claims the benefit of Provisional 
Application No. 60/116,362 filed January 19, 1999, 
incorporated herein by reference. 

BACKGROUND OF THE INVENTION 

This invention is concerned with signal processing of 
transients, and more particularly, with multiexponential 
signal processing. 

In the recent REVIEW article entitled "Exponential 
analysis in physical phenomena". Review of Scientific 
Instruments, Vol. 70, No. 2, Feb. 1999, pages 1233-1257 
(incorporated herein by reference) , authors Andrei A. 
Istratov and Oleg F. Vyvenko consider various aspects of 
multiea^onential analysis. It is apparent from the text of 
the REVIEW article, and from the over 300 literature 
citations, that exponential amalyeis in physical phenomena 
is a subject of considerable interest in many scientific and 
t cbnological disciplines, including, inter alia, solid 
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state physics, medicine, biology and biophysics, geophysics, 
optics, engineering, chemistry and electro-chemistry. 

As pointed out in the REVIEW article, many physical 
phenomena are described by first-order differential 
equations whose solution is an exponential decay. The 
amplitude and time constant of the exponential decay carry 
information about the nature of the phenomenon being 
studied ♦ As commonly happens, a number of exponential 
processes take place simultaneously, and equipment employed 
in analyzing such exponential processes (multiexponential 
decays) yields a signal which is a sum of a plurality of 
exponential components • 

Obtaining useful information from multiexponential 
decays is not a simple task. Many methods of 
multiexponential analysis and many algorithms for this 
purpose have been proposed, but all of them have 
limitations • 

U.S. Patent No. 5,517,115 issued May 14, 19S6 to 
Manfred G. Prammer (incorporated herein by reference) 
proposes a method and apparatus for efficient processing of 
nuclear magnetic resonamce (HMR) echo trains in well 
logging. A priori information about the nature of the 
expected signals is used in an attempt to obtain an 
approximation of a model using a set of pre- selected basis 
functions. A singular value decomposition (SVD) is applied 
to a matrix incorporating information about the basis 
fxinctions, and is stored off-line in a memory. During 
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actual measurement, the apparatus estimates a parameter 
related to the signal-to-noise ratio (SNR) of the received 
NMR echo trains and uses it to determine a signal 
approximation model in conjunction with the SVD of the basis 
function matrix. This approximation is used to determine, 
in real time, attributes of the earth formation being 
investigated. 

Techniques such as those described in the Prammer 
patent rely heavily upon a priori information about the 
nature of the expected signals. In attempting to obtain a 
reliable estimate of an unknown model, i.e., parajcieters of a 
subject being investigated, such techniques rely on one or 
more numerical algorithms for multi exponential analysis 
which use so-called "fitting routines" in an attempt to fit 
data to the model. 

Underlying the present invention is the discovery that 
substantially inqproved results can be attained withotxt using 
fitting routines, and, instead, by emphasizing better linear 
resolution. As pointed in the REVIEW article, fitting 
routines work well only if the hypothesis of the number of 
multiexponential components is correct and an initial 
approxixaization is close to the true solution. The present 
invention, which does not use fitting routines, avoids the 
problems that are inherent in the use of such routines. 



wo 00/43906 PCT/raoO/00212 



BRIEF DESCRIPTION OF THE INVENTION 

An object of the present invention is to provide 
improved methods and apparatus for the analysis of 
transients and for obtaining useful information therefrom, 
5 More particularly, an object of the invention is to 

provide improved multi exponential signal processing. 

One aspect of the invention involves a computer- 
readable medium containing a set of coefficients that define 
a transform operator such as a matrix. 
10 Another aspect of the invention involves a method of 

calculating a transform operator utilizing a plurality of 
resolution functions. 

Another aspect of the invention involves a method of 
multiexponential signal processing in which multi exponential 
15 signals are san^led, and the above-mentioned transform 
operator is applied to the ScLS^led signals. 

Another aspect of the invention involves an apparatus 
for multiexponential signal processing that coxi^rises a 
signal processor including the above-mentioned transform 
20 operator. 

The present invention starts with the const^ruction of 
appropriate transfozra operators. Once appropriate transform 
operators have been constructed, they are incorporated in 
signal processors of analytical instrumentation for 
25 processing data. Generally, such instrumentation includes a 
computer, as is well known in multiexponential analysis* 
Multiexponential data signals are sensed or detected by 
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conventional equipment and are input to the transform 
operator of the computer for signal processing. The signals 
may be applied in real time or they may be read out from a 
suitable storage medium* 
5 Typically, the signals are applied to the transform 

operator in digital form after conventional sampling and 
analog- to-digital conversion. For example, digital samples 
of multiexponential decays may be obtained at eq[ually- spaced 
instants in time, beginning at or just afte^ the start of a 

10 multiexponential signal. 

The invention makes use of linear resolution to obtain 
a better estimate of an iinknown model. The present invention 
provides the very desirable property of optimal linear 
resolution of the unknown model when plotted against the log 

15 of the time constant of the decay curves. 
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BRIEF DESCRIPTION OF THE DRAWINGS 

The invention will be further described in conjunction 
with the accompanying drawings, which illustrate preferred 
(best mode) emboddLments , and wherein: 
5 Fig. 1 is a diagram showing a matrix with coefficients 

^11 * * * ^or transforming data d,^ • • • d^ to produce 

parameters m^^ . . . m,^ of an estimate of an unknown model; 

Fig. 2 is a flow chart showing a method for calculating 
the coefficients of a row of a transform matrix in 
10 accordance with the invention; 

Fig. 3 is a table showing, inter alia > coefficients for 
three rows of a matrix, for r's of interest, namely 1, 10, 
100; 

Fig. 4 is a diagram showing data fxinctions, resolution 
15 functions, noise response, ajxd point spread functions for a 
matrix constructed in accordance with the invention where 
the noise gain (NG) is 1.0000; 

Fig* 5 is a similar diagram for another matrix 
constructed in accordance with the invention for a noise 
2 0 gain of 3.1623; 

Fig. 6 is a diagram showing resolution ftinctions for 
four matrices constructed in accordsuice with the invention 
with different noise gains; 

Fig. 7 is a diagrcun showing point spread functions 
25 associated with four matrices constructed in accordance with 
the invention for different noise gains; 
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Fig. 8 is a block diagrsun showing apparatus for 
processing data in accordance with the invention; 

Fig. 9 is a diagram showing decay curves of an MRI; and 
Fig. 10 is a diagram showing MRI relaxation 
distributions for four matrices constructed in accordance 
with the invention with different noise gains. 

DETAILED DESCRIPTION OF THE INVENTION 

One of the principal objectives of the present 
invention is to provide signal processing of transients, 
such as multi exponential decays, producing outputs that are 
better estimates of an \inknown model to be investigated. 
Such estimates permit better interpretation of data, so that 
a user (researcher, physician, scientist, or engineer, for 
example) can obtain more accurate information as to the 
nature of the unknown model. Considered from one point of 
view, the invention may be looked upon as a better digital 
lens that provides improved resolution, just as a better 
optical lens provides improved resolution. 

A first step in achieving the objectives of the 
invention is to construct an improved trsuisform .operator, 
conveniently in the form of a matrix. The mcuiner in which a 
transform matrix may be constructed, pursuant to the 
invention, will be described in detail later. The actual 
transform operator will depend upon its intended 
application, for example, medical imaging or well logging. 
For any application, several different transform operators 
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may b constructed, to provide a user with greater 
flexibility. 

The present invention requires an understanding of what 
is referred to in the art as estimating a solution of a 
linear inverse problem. The linear inverse problem is one 
of communicating to an interpreter what is known and what is 
not known about an xmknown model. The almost universal 
practice in the prior art for estimating the solution of a 
linear inverse problem is to calculate one or more of an 
infinite number of estimates of the model which fit the 
data, i.e., which reproduce the data to within the noise 
that is present. Whenever possible, a priori information is 
used to choose which of the estimates to calculate. 
However, a priori information may not be sufficiently 
available or may be suspect. 

In applying the present invention to multiexponential 
decays, an estimate of an unknown relaxation distribution 
(model) is obtained by linearly resolving each point of the 
xmknown relaxation distribution as precisely as possible 
within the limits of the noise • In general, such estimates 
do not reproduce the data. Rather, they are obtained by 
optimizing the linear resolution in a manner that will be 
described later. 

The term "linear" used herein in connection with 
"resolution" has the following meaning: 
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Consider a transform operator A which, transforms 
functions (or vectors) xl and x2 to yl and y2 , respectively. 
As equations, this is stated: 

A xl = yl and A x2 = y2 • 
The traxisform operator A is defined as "linear" if for any 
two real nviznbers a and b 

A (a xl + b x2) = a yl + b y2 . 

A transform operator which is a matrix will always have 
linear resolution. It is possible to construct non-matrix 
transform operators which behave similarly to matrices for 
only a restricted set of functions (or vectors) . Moreover, 
a matrix can be derived from a non-matrix transform operator 
that expresses this behavior. 

The present invention involves the following 
relationship between a true model m'(y) and data dj^ (e.g., 
multi exponential decays) : 

where I is the appropriate domain of definition and 
g,^(y) represents a data function, more particularly, one of 
the N data functions which coxc^ose the data kernel of a 
linear transform. A data function is a mathematical 
representation of how a model is mapped to a particular data 
point. Equation (1) is referred to as "the forward 
problem. " 
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Decaying systems generally have a point, usually 
defined as t=0 (where t is time) , at which a system begins 
decaying. The system decays to a constant value, o£ten 
zero, as t approaches infinity. In nuclear magnetic 
resonance (NMR) , for example, the t=0 point is when an 
excitation pulse energizes a sample. 

The multiexponential forward problem commonly used in 
data analysis and inverse theory has the form 

In this equation di^ represents the data; m^ represents the 
unknown model, and -e^Wn represents the data function, 
where r designates the time constant of the exponential 
decay. 

Equation (2) can be expressed in the continuous form 

dk = f^T^^^ir - n) e-^f^ (3) 

where 6() is the Dirac delta function. 

Incidentally, as used herein, the range of indices such 
as i, j and k are determined by the equation in which they 
are used. For example in the equation 
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the index i wou^d be assumed to remge over all the data 
functions ^(y) • which are normally numbered with positive 
integers starting at 1. But in the equation 

NGi^^^j (5) 

the index i ranges over the rows of a transform matrix and 
the index j ranges over the columns of the transform matrix, 
which is the same range as the data functions in the 
previous equation. 

It is desirable to plot the relaxation distribution 
(unknown model) versus the logarithm of r, and to determine 
the resolution of a relaxation distribution on a log scale. 
Applying a change of variables y=ln(T> to equation (3), 
without sizr^lif ication, yields 

dk = f^lZ^H^ - e^Oe-***'^(c»'dy) (6) 
Applying the Dirac delta function identity 

= IZjJ^^i^ " ^) each f(xi) =0 (7) 
to equation (6) along with stamdard simplification yields 

Substituting 

i 
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into qua t ion (8) gives 



(10) 



where m^(y) represents the xinknowa model, and c"'^**^ 
represents the data function, which may be expressed 

as : 



Hereinafter, the form of the multiexponential transform in 
equation (10) will be termed the forward problem* 

It should be noted that the data function approaches 1 
as T approaches infinity for all values of t,^* Pursuant to 
the invention^ the data function of equation (11) has been 
found to be very useful in multiexponential analysis, but 
other data functions may be appropriate for other 
applications of the invention. 

As stated earlier, a first step in achieving the 
objectives of the invention is to construct a transform 
operator, such as a transform matrix, which maps data (e.g. 
decay signals) to xinknown model parameters. Fig. 1 Is a 
diagram showing a matrix with coefficients a^^^ . , . a^^. for 
transf orzoing data d^. . • d„ to produce parcuneters m^ • • • xa^ 
of an unknown model. The transform matrix is typically a 
matrix in which each row corresponds to a r o£ interest. 
Selection of r's of interest and data points for initial 
coefficients will be guided by available inforzaation in the 
particular field in which the matrix is to be used. 



(11) 
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In accordance with the invention, it is desired to 
obtain an estimate o£ an unknown model with linear 
resolution, and since the data fiinctions used are linear 
functions of the model, the estimate of the unknown model 
can be calculated by multiplying the data by a matrix. The ^ 
matrix is chosen so that each point of the estimate of the 
unknown model linearly resolves the corresponding point of 
the model as well as possible with an acceptable noise gain 
i,e., with optimal linear resolution. Since matrix 
multiplication is a linear operation, it yields an estimate 
which does not necessarily reproduce the data, but which 
does have linear resolution. Linear resolution is a 
desirable property of an estimate, because each point of the 
estimate resolves the corresponding point of the model in 
the same way, independent of any particular model. 

A goal in constructing a trajisform matrix in accordance 
with the invention is to calculate linear combinations of 
the data fxinctions which yield a resolution function that 
resolves as small a region of the uziknown model as possible. 
Accordingly, it is preferred that each resolution function 
corresponding to a row of the matrix be characterized by 
optimal linear resolution. Preferred criteria for selecting 
coefficients of a transform matrix which yield optimal 
linear resolution are described later. Together, the 
resolution function and its noise gain give a concise 
formulation of the ambiguity of a point in an estimate of an 
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unknown model from information given by data and 



corresponding data functions. 



Calculation of a transform matrix proceeds row by row. 



It is convenient to use tbe same noise gain for all rows of 
a transform matrix, but this is not required. As noted 
above, each row normally corresponds to a t of interest. 
The rows can be ordered by increasing r from top to bottom 
of the matrix. A spacing of 16 r values per decade has been 
found to work well. As an example, r values may be 
calculated between O.lt^^and lOt^, where t^„ is the 
smallest time at which a decay signal is to be sampled and 
^nax largest time. Each row of the matrix corresponds 

to a resolution function that is centered on a r of 
interest. 

It is presently believed that the best way to calculate 
the coefficients of a particular transform matrix (which, 
incidentally, may have only one row) is to solve a 
constrained minimization problem. 

The information needed before the constrained 
minimization can be performed are (1) the data function, 
gi(y>* for each data point, d^, (2) the stcuidard deviation of 
the expected noise, a^, of each data point, (3) the desired 
noise gain, NG, and (4) the r of interest, t^, on which the 
resolution function is to be centered. 

The constrained minimization problem to be solved is to 
minimize 1 in the equation: 




(12) 
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by varying b^, where 

<4) 

i 

and where b^ are trial coefficients of a transform matrix 
row and R(y) is a trial resolution fiinction. More 
information on resolution functions can be found in Parker, 
Robert L. 1994; Geophysical Inverse Theory ; Princeton, New 
Jersey; Princeton University Press, incorporated herein by 
reference* 

The trial resolution function preferably complies 
substantially with the following constraints: 

(i) R(y) > 0 [Nonnegative constraint] 

(ii) R(y,,) > 1 [Peak constraint] 

(iii) R(y)monotonically decreases from [Monotonicity 
constraint] 

^^^^ . £6<^3 ^ l^bigiiy)dy (Noise gain constraint] 

Constraint (iii) may be satisfied automatically for the 
data fxmction used for the multiexponential problem- Thus, 
while it is stated as a separate constraint, it is to be 
xmderstood that it may, in fact, be satisfied inherently. 

To arrive at the final coefficients, a^^, for row i of 
the transform matrix, once the trial coefficients b^ which 
minimize I have been found, it is highly desirable that the 
trial coe££i.ci.ezit8 be scaXed by a constEUXt so tbat tlie airea 

15 
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of their associated resolution function is unity on a log 
scale. This is accomplished by 

= V bk9kiy)dy ( 13 ) 

Making the area of the resolution function unity insures 
that each point in the estimate of the unknown model is a 
local average. 

Incidentally, for a data point d^,, sampled on a decay 
curve at time t^, the data function previously described' can 
be modified by a balancing function B (y) , so that the data 
function becomes 

This modification of the data function has the effect of 
making the resulting trcuasform matrix generate an estimate 
of the unknovm model multiplied by B (y) but with the 
recjuired noise gain. While many balancing functions are 
possible, a useful balancing fxmction ie = , The 

choice w s 1 may be useful since it increases the amplitudes 
of the relaxation components at late time constants by a 
factor of T. 

Constructing a transform matrix in which each 
resolution fiinction corresponding to a row of the matrix is 
characterized by optimal linear resolution preferably 
involves a process termed constrained optimization. In the 
preferred embodiment of the invention, this process includes 
an outer loop, a middle loop and an inner loop. The outer 

16 
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loop converts the constrained optimization problem to a 
series o£ iincons trained multidimensional minimization 
problems using the simple penalty method described by 
Fletcher (Fletcher, R. 1987 Practical Methods of 
5 Optimization ; Toronto; John Wiley & Sons) incorporated 
herein by reference. The middle loop converts the 
tincons trained multidimensional minimization problems into a 
series of one dimensional (ID) minimization problems via the 
conjugate gradient method described by Press et al. (Press, 

10 W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P. 
1992; Numerical Recipes in C; The art of scientific 
computing ; 2nd Ed,; Cambridge UK; Cambridge University 
Press) incorporated herein by reference. The 1 dimensional 
minimization problems are solved using the golden section 

15 search or parabolic interpolation described by Press et al. 

A desirable condition for each resolution function is 
that its area should be svibstantially minimized. 
Calculating the area of the trial resolution function X to 
be minimized, given the trial coefficients, b^, can be 

20 accozE^lished using the equation: 

I^i:biAi (15) 

t 

where is the area of the data function, gi (y) / given 
by the equation 
25 . rn 



^^f Si(y)dy 
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The integral's upper limit of infinity has been replaced by 
the value Yl, which is a large finite ntamber. Experience 
has shown that Yl = ln(1000*t^), is a good choice. 

The simple penalty method converts the constrained 
optimization problem to an unconstrained one via residues, 
r„. A residue is a measure of how much a particular 
constraint is violated. If a particular residue is greater 
than zero, then the corresponding constraint is considered 
satisfied. The more negative the residue, the more the 
constraint is violated. The constrained optimization 
problem then becomes the unconstrained one of minimizing I** 
in the equation: 

i m 

where min(x,y) returns the minimxim and x and y* The 
value of P starts small, and is then increased by a small 
factor, for example, 0.1, after each unconstrained 
minimization. The solution to one constrained minimization 
is used as a starting point for the next constrained 
minimization. The value of P is increased until I** has 
stabilized. The test for stability is whether the value of 
I** changes by an accumulative factor of 10"* over four 
consecutive increases in P. 

Each of the non-negative, peak, monotonicity and noise 
gain constraints will usually have associated residues. 
Before the residues Ccoi be calculated, the data functions 
and the resolution functions are discretized. 

18 
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The data functions and the resolution fimctions are 
discretized by sampling them at e.g., 16 points per decade 
of T starting at t= 0.01 time units and continuing to r= 
1000*t^ time units. The discretized data functions, gi (y) , 
are represented as g^j . The total niunber of points at which 
the functions are sampled is denoted L. The discretized 
trial resolution function is represented by . 

Thus 

Ri = Tl^iSu (18) 

i 

The residues of the noiuiegative constraint, r^*^, are 
defined to be 

rf^'-^Mt (19) 
giving £< residues. 

The residues o£ tb.e peak constraint « r^, are defined to 

be 

r^^Rt^-l (20) 

where Ip is the index of the discretized trial 
resolution function at the r of interest. This yields one 
residue. 

The residues of the mono tonicity constraint are defined 

to be 

r/^^rrKi-A-i 2>l>lp (21) 

and 

19 
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L-1 



(22) 



These two equations together yield another L - 1 residues. 

The residue of the noise gain constraint is defined to 

be 



This equation yields one additional residue. 

In total there are 2*L+1 residues to express the four 
constraints , 

The conjugate gradient technique is explained by Press 
et al. The termination conditions used may be an 
acciimulative factor change in I*" of no more thcin 10"'' over 2 5 
consecutive iterations. The total number of iterations can 
be limited to 3000*N, where N is the number of data points. 

The one dimensional m±ndLmization problems are solved 
using the golden section search or parabolic interpolation 
described by Press et al. The ID optimization is terminated 
If the functional value does not decrease by a certain 
fraction each step. The fraction can be set to 10'^, for 
example. A step limit of 1000 on each ID optimization can 
be set to prevent infinite or unproductive near infinite 
loops • 

The linear transform for the trial row coefficients to 
the trial resolution function would usually be computed a 
large nximber of times during each one dimension optimization 
and would be computationally demanding. However, since the 
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relationship between the trial row coefficients and the 
trial resolution function is always linear, this linearity 
can be used to greatly reduce the nvunber of times the trial 
resolution functions have to be calculated directly from the 
trial coefficients . 

In one dimensional optimization, a scalar parameter c 
is varied along a line with direction d^; the direction 
being provided by the conjugate gradient method. Therefore, 
for the trial coefficients b^, optimization occurs along a 
line defined by 

bn=bf+cdt. (24) 

where b^ is also supplied by the conjugate gradient method* 
Because of the linearity relationship between b^^ and R^, a 
similar equation can be written for the corresponding trial 
resolution fxinction R^, 

Ri=:Ji}+cd^ (25) 

where 

^~J2^9u (26) 

i 

and 

^^^J2^9u. (27) 

Therefore, along each line of optimization, the linear 
transform from the trial coefficients to the trial 
resolution function only needs to be calculated twxc , to 
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calculat:e and d^, instead of once for every value of c 
considered. This results in a many fold reduction in the 
one dimensional optimization time. 

The effectiveness of the minimization can be checked by 
how close the value of the peak of the resolution function 
is to unity, since the minimxim area should yield a peak 
value of exactly unity. 

An upper limit on the noise gain, NG, of a particular 
point in an estimate can be applied by the noise gain 
constraint stated earlier, where the area of the resolution 
function is included because the area of the resolution 
function must be normalized (set equal to 1) before the 
noise gain is calculated. 

Figure 2 is a flow chart for calculating the 
coefficients of a row of a transform matrix pursuant to the 
foregoing description, using a digital computer 
conventionally. Coefficients for successive rows of the 
matrix can be calculated in the same manner for each r of 
interest and for each noise gain deemed appropriate. For 
calculating the coefficients of any row, the coefficients of 
the preceding row can be used as starting points*.. 

Constrained optimization methods may perform better if 
all the data have noise with a standard deviation of 1. 
Fortunately, it is quite sample to modify the data functions 
so this is the case. Given the real data functions, gi (y) # 
the adjusted data function, gi*(y) is 

fl^(y)=^<(v)M (28) 
22 
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The adjusted data functions are then supplied to the 
constrained optimization method along with the assumption 
that the standard deviations of noise at all the data points 
are 1. The final coefficients, a^^, produced by the 
constrained optimization based on the adjusted data 
functions, need to be corrected for the adjustment. The 
correction to the final coefficients is 

= (29) 

If a data set has 100' s, 1000' s, or even more evenly 
spaced points on a decay curve, it is much more efficient 
computationally for calculating the coefficients of the 
transform matrix (and also for applying it to data) to 
average adjacent data points together to create a new group 
data point. The corresponding data functions must also be 
averaged together to get the grouped data function 
corresponding to the new grouped data points. The size of 
the groups is importouit. While, in general, it is best to 
have larger groups at later sample times, groups which are 
too large will reduce the resolution of the resolution 
functions. Groups which are too small are inefficient. A 
logarithmic group size appears to be a good choice. 
The equation 

Gn = max(l, int{C 10^^^). (30) 

appears to produce good group sizes for C = 0.1, D = 10, 
where 1 is a positive Integer starting at 1 and increasing 
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to whatever value is appropriate • The function xaax(x,y) 
returns the maximum of x and y, and the function int(z) 
rounds 2 down to the nearest integer. For 512 data points 
the above equation gives group sizes of 

5 I 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 3 3 5 6 7 10 12 15 19 25 31 39 50 63 79 100 30 

Occasionally/ instead of measuring data at a point in 
time on a decay curve, data will be measured by averaging 
over a window between two points in time. In this case, the 
data function can be approximated by averaging together a 

10 larger number of sample points over the window. 

Modification of the data functions will also be 
required if a trigger that initiates a decay curve is not a 
single point in time. For exan^le, in time resolved 
fluorescence spectroscopy, the flash of light triggering the 

15 decay will last a finite length of time. If the intensity 
versus time for the flash is Ii(t) then each data point will 
be convolved with this function. The corresponding data 
functions will also have to be convolved with the same 
function* 

2 0 Fig. 3 illustrates a trcoisform matrix constructed in 

accordance with the invention. Three sets of coief f icients 
for T=l, 10 and 100 are shown as columns so that they fit on 
the page. Each set has 48 coefficients calculated from 48 
data functions. The first 32 data fxinctions correspond to 

25 32 saxi^ling points that are evenly spaced by one time unit 
at times 1 to 32. The next 16 data functions are evenly 
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spaced by 30 tdLme units between time 62 and time 512. The 
standard deviation of the noise for the first 32 points is 
assumed to be 1, and the standard deviation of the noise of 
the next 16 is assumed to be 0.18257. The standard 
5 deviation of the noise of the last 16 is a factor of 
l/sqrt{30) less than the first 32. This drop in the 
standard deviation of the noise would result if the cut-off 
frequency of a low pass filter before the analog- to-digital 
converter were dropped by a factor of 3.0 before the point 

10 was acquired. For the sake of simplicity, the matrix was 
constructed without the use of balancing functions or data 
function groups. 

Fig. 4 illustrates data fxinctions, resolution 
functions, noise response, a^d point spread functions for a 

15 matrix constructed in accordance with the invention where 
the noise gain is 1.000. The data function curves are for 
time samples at t=l, t=:2 . . . t=N from left to right. Each 
resolution fiinction corresponds to a set of data functions. 
The abscissa in each diagram is in units of t on a log scale 

20 [ln(T)], and each resolution f\inction is localized about a 

particular r value. As a result, each point spread, function 
tends to be localized aO^out a particular r value. 

Fig. 5 is a diagreun similar to Fig. 4 but for a matrix 
constructed with a noise gain of 3*1623. It should be noted 

25 that in each of Figs. 3, 4, and 5 there are 48 data points 
(time samples) . 
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A comparison of resolution functions produced by 
inatrices with different noise gains is shown in Pig, 6. 
Higher resolution is achieved with higher noise gains, but 
there is a trade-off between resolution and noise gain. 
Greater noise tends to make the results achieved less 
reliable . 

In addition to resolution functions, performance of a 
transform matrix can be judged using the PSF's and noise 
response. To judge the performance of a transform matrix, 
information is required as to corresponding time points on 
the decay curve at which data should be acquired as well as 
assumed noise at each data point. 

Fig. 7 shows point spread functions associated with 
four matrices constructed in accordance with the invention 
for different noise gains. PSF's can be calculated using 
simulated data at required time points for monoexponential 
decays at t = 1,2,5,10,20,50,100, ... time units, assuming 
the initial amplitude of the decay curve is 1. No noise is 
added to the simulated data. It should be noted that in 
each of Figs. 6 and 7, and also on Figs. 9 and 10 to be 
described, there are only 32 equally spaced data points 
(with the same noise) . 

To calculate the noise response, several realizations 
of a decay curve which are purely noise are generated. The 
noise may be assumed to be Gaussian. The standard deviation 
of the generated noise is scaled so that the standard 
deviation of the noise of the first data point is 1. Then 
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the sampled noise decay curves are multiplied by tlie 
transform matrix to obtain the relaxation distribution. 
Several realizations of the relaxation distribution of the 
noise can be plotted to obtain a "feel" for the 
distribution. 

An important characteristic at each point in the 
estimate of the unknown model is the standard deviation due 
to the noise in the data. If the noise in the data is 
iincorrelated/ has a mean of zero, and all points have the 
same finite standard deviation, it can then be characterized 
by a single standard deviation. The noise gain for each 
point can be defined to be the standard deviation of the 
point divided by the standard deviation of the noise in the 
data and can be calculated directly from the coefficients 



Once a transform matrix has been constructed, 
experiments can be performed in which the spacing of the 
data acquisition points, the number of data acquisition 
points, £uid the grouping of data acquisition points are 
changed. The results of these experiments can be compared 
by observing the resolution fimctions and/or the PSFs that 
are produced. As a rule of thximb, the linear resolution is 
proportional to the height of the resolution function 
provided that the resolution function has \mity area on a 
log scale. 




(31) 
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Each trial resolution function is localized about a 
point of interest in the vmknown model by requiring the peak 
of the resolution function to have a value greater than 
unity at a T of interest and then by minimizing the area of 
5 the resolution function. The effectiveness of the 

minimization can be checked by how close the value of the 
peak of the resolution function is to iinity, since the 
minimum area should yield a peak value of unity. 

From the foregoing description, it is evident that each 

10 transform matrix is constructed so that each point of the 
estimate of the xinknown model linearly resolves the 
corresponding point of the unknovTn model as well as possible 
within an acceptable noise gain. Since matrix 
multiplication is a linear operation, it yields an estimate 

15 which does not necessarily reproduce the data but does have 
linear resolution. With linear resolution, each point of 
the estimate resolves the corresponding point of the xinknown 
model in the same way independent of any particular unknown 
model. A highly desirable property of linear resolution in 

20 providing an estiioate of an unknown model is that it enables 
a human interpreter to obtain an intuitive "feel" of what 
the data reveals and does not reveal alDOUt the unknown 
model . 

Each resolution fvmction gives a concise mathematical 
25 description of the linear resolution at each point of the 
unknown model and is independent of the unknown model and 
the data. Viewing the treaisform matrix as a digital lens. 
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linear combinations o£ the data functions are calculated 
wh.ich yield a resolution fiinction that resolves as small a 
region of the unknown model as possible. 

After a transform matrix is constructed and selected, 
5 it is incorporated in a signal processor of a computer, as 
software or hardware, for example, as indicated in Fig. 8. 
In this figure the INPUT represents a source of sampled 
digital signals such as multiexponential decays, e.g., NMR 
decays obtained from well - logging . The DATA PROCESSOR and 

10 STORAGE are components of a conventional digital computer. 
The OUTPUT MODULE may have conventional DISPLAY, NETWORK 
INTERFACE, and PRINTER COMPONENTS, for example. 

An estimate of an unknown model is calculated by 
multiplying data, e.g., multiexponential decays, by the 

15 transform matrix. Several transform matrices for different 
noise gains, for example, can be incorporated in the signal 
processor and accessed selectively to provide a user with 
greater flexibility. 

As stated earlier, an object of the present invention 

20 is to provided a better estimate of an \mknown model, from 
which useful information can be obtained. In the 
aforementioned Prammer patent. Figure 8 of the patent is a 
mapping of estimated NMR signal decay times into pore sizes 
of an investigated earth formation. The curve of Fig. 8 is 

25 an estimated relaxation distribution. The present invention 
provides a better estimate of the relaxation distribution, 
and also a better signal-to-noise ratio (SNR) . Resolution 
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functions produced by matrices constructed in accordance 
with the invention are superior, in terms of peak value and 
localization to resolution functions produced by matrices of 
the Prammer patent for the same data and noise. 

In general, the present invention can be used to 
provide better estimates of an unknown model than the prior 
art, estimates from which more relicd^le and useful 
information can be obtained. These improved results are 
achieved by emphasizing linear resolution irrespective of 
whether data fits a model, and, in fact, without any attempt 
to fit data to a model. The invention is particularly 
useful in multiexponential signal processing, such as in the 
analysis of multiexponential decays* 

As mentioned earlier, witli particular reference to the 
REVIEW article, multiexponential analysis is useful in a 
host of scientific and technological applications. One of 
those applications, NMR, such as NMR in well -logging, has 
already been discussed. Another application is MRI 
(magnetic resonance imaging) . Fig. 9 shows decay curves in 
an. MRI application. The decay curves represent pixels taken 
from a series of 32 magnetic resonance images of the brain 
of a multiple sclerosis patient. The Tj relaxation data 
were acquired with a 10ms sample time out to 320ms. The 
strength and decay of the signals contain valuable 
information about the tissues. Some of the major tissue 
types of interest in multiple sclerosis are normal appearing 
white matter (NAWM) , c rebral spinal fluid (CSF) and 
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lesions, which can be classifi d as chronic or acute. Fig. 
10 shows relaxation distributions yielded by applying to the 
decay curves, of Fig. 9 transform matrices of the invention 
with various noise gains. The larger the noise gain of the 
5 estimate, the wider the range of relaxation estimates 

available. This is a characteristic which shows up in the 
resolution functions and is due to the fact that the larger 
the noise gain the more likely there is a feasible 
resolution function available . 

10 Estimating the noise in the estimates shown in Fig. 10 

can be accomplished in several ways. The first is to 
estimate the noise in the data and multiply by the noise 
gain for each transform matrix. The ideal way to measure 
the noise in the data is to repeat the measurements a large 

15 number of times and calculate mean, standard deviation euxd 

covariance. These statistics can then be propagated through 
transform matrices using standard statistical procedures. 
Unfortunately, the measurement of the decay curves takes 
about 20 minutes to con^lete on a patient, so large mmbers 

20 of repetitions are impractical. 

Another way to measure the noise in the data i@ to 
estimate the standard deviation from previous measurements 
using the same instrximent. This is not always reliable 
because the noise can vary from sample to sample. For the 

25 data in Fig. 9, previous measurements of noise gives the 
standard deviation of about 10 scanner units. After 
multiplying the standard deviation by the noise gain, the 
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predicted standard deviation for the noise in Fig. 10 is 
100. 

A third way to estimate the noise is to consider that 
while the noise gain increases by a factor of 10 in Fig. 10, 
the resolution gain increases by only 1.4 for relaxation 
rates around five sample times. In Fig. 10, a portion of 
the signal between 3 0 and 2 0 0ms increases proportionately to 
the noise gain. This strongly suggests that the portion is 
due to random uncorrelated additive noises in the data. It 
is possible then to measure the standard deviation of the 
noise between 3 0 and 2 0 0ms and work back to the noise in the 
data. 

Further applications of the invention will be apparent 
to persons of ordinary skill in the art. For example, 
decaying sinusoids are a common problem in inverse theory. 
A decaying sinusoid can be hajadled if it is band pass 
filtered at the particular bandwidth of interest and then 
the magnitude of the decay curve taken. The relaxation 
distribution of the magnitude decay curve caji. then be 
calculated. 

Any row of a transform matrix, since it corresponds to 
a resolution function, can be applied to a time series in 
the same way as digital filters. The "relaxation" digital 
filter will be useful for applications such as ultrasound 
and radar. Using qn^adrature detection, it would be possible 
to measure the magnitude of a reflection off o£ an 
interface. If the signal from an interface of interest 
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oscillated for a whil after th initial sound wave had 
passed, the decay time of the oscillation would reveal 
information about the interface. Applying various 
relaxation digital filters to the ultrasound time series 
5 would allow characterization of the interfaces. If it were 
desirable to detect a particular range of relaxation time, 
selection of the coefficients of a relaxation digital filter 
would give a resolution function with a more boxcar- like 
shape . 

10 Other applications of the invention include 

photoluminescence and time -resolved fluorescence 
spectroscopy of biological and other types of samples, 
electrical signals radiating from ore bodies in geophysical 
exploration and acoustic, electrical and electromagnetic 

15 decay processes. Xt is possible to design a transform 
matrix which handles data integrated over a window or 
unevenly spaced in time by modifying the data functions. 

Principles of the invention can be applied to problems 
which data functions other than decay curves. If the data 

20 functions are cosine ftmctions, resolution functions can be 
generated for low pass, band pass and high pass filters as 
well as windows for discrete Fourier transforms. A low pass 
filter ccm. be designed by requiring maximum area below the 
cutoff frequency, minimum area above the cutoff and an 

25 optional req[uirement of monotonicity to eliminate wiggles. 

The bounds on all parts of the filter would be 0 and 1. A 
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limit on broadband noise gain could also be imposed to 
improve tbe robustness of the filter. 

While preferred forms of the invention have been shown 
and described/ these forms are intended to be exemplary, not 
restrictive, and it will be apparent to those skilled in the 
art that modifications can be made without departing from 
the principles and spirit of the invention, the scope of 
which is set forth in the appended claims. 
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WHAT IS CLAIMED IS: 

1. A fcomputer- readable medium containing a treuasfona 
operator constructed to provide an estimate of an unknown 
model with substantially optimal linear resolution. 

2. A computer- readable medium according to Claim 1, 
wherein the transform operator comprises a matrix having at 
least one row of coefficients corresponding to a resolution 
function* 

3. A computer- readable medium according to Claim 2, 
wherein the resolution function has substantially no 
negative values • 

4. A computer-readable medium according to Claim 2, 
wherein the resolution function has a peak value of at least 
substantially 1 substantially at the point of interest. 

5. A computer -readable medium according to Claim 2, 
wherein the resolution function monotonically decreases from 
its peak value. 

6. A coH^uter- readable medixm according to Claim 2, 
wherein the resolution function has sxibstantially no 
negative values, has a peak value of at least substantially 
1 substantially at the point of interest, and decreases 
substantially monotonically from its peak value. 



wo 00/43906 



PCT/IBOO/00212 



1 7. A computer- readable medium according to Claim 6, 

2 wherein the resolution function complies suibs tantially with 

3 the following noise gain constraint: 

"if ^ i^Hw ^ 

\ 53 ^^^^4 Yl^i9iiy)dy [Noise gain constraint] 

1 8, A computer-readable medixim according to Claim 2, 

2 wherein the matrix has a plurality of rows of coefficients, 

3 each row corresponding to a different resolution function. 

1 9. A computer-readable medi\un according to Claim 8, 

2 wherein the matrix is constructed for use in 

3 multi exponential decay signal processing and wherein each 

4 resolution function is substantially centered on a different 

5 time constant. 

1 10. A computer- readable medium according to Claim 9, 

2 wherein each coefficient has a corresponding data fiinction 

3 in the form c"***"*' or derivations of this form, where y is 

4 the natixral logarithm of a time constant and t„ is a 

5 saxapling point on the decay signal. 

1 11. A coii5)Uter- readable medium according to Claim 2, 

2 wherein the coefficients are selected so that the 

3 corresponding resolution function has a peak value of at 

4 least substantially 1 substantially at the point of interest 

5 and wherein the area of the resolution function is 

6 substantially minimized. 
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12. A computer -readable medixun according to Claim 11, 
wherein the area of the resolution fiinction is set to be 
substantially 1 on a logarithmic scale* 

13. A method of constructing a transform operator in 
which a plurality of coefficients are selected to produce a 
corresponding resolution function that provides 
substantially optdLmal linear resolution. 

14. A method according to Claim 13, wherein the 
resolution function has substantially no negative values, 
has a peak value of at least substantially 1 siibstantially 
at the point of interest, decreases substan.tially 
monotonically from the peak value, and has a gain within a 
predetermined range* 

15. A method according to Claim 13, wherein the area 
of the resolution function is set to be substantially 1 on a 
logarithmic scale. 

16. A method according to Claim 13, wherein the 
transform operator is constructed to provide substantially 
optimal linear resolution between data and an unknown model 
without any consideration of whether the data fits the 
xinknown model. 
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1 17 . A method o£ znultiexpon ntial signal processing, 

2 which comprises: 

3 sampling a multiexponential signal and applying the 

4 transform operator of Claim 1 to the sampled signal. 

1 18. Apparatus for multiexponential signal processing, 

2 which comprises a signal processor that has the transform 

3 operator of Claim 1 . 

1 19 . A method of constructing a transform operator, in 

2 which a plurality of coefficients are calculated using data 

3 fxinctions that emphasize linear resolution irrespective of 

4 data . 

1 2 0. A method according to Claim 19, wherein each data 

2 function is in the form c"***"*' or derivations of this form, 

3 where y is the natural logarithm of a time constant and tj, 

4 is a sampling point on a decay signal. 

1 21. A method of exponential signal processing, which 

2 comprises: 

3 providing a sampled multiexponential signal; and 

4 applying the scunpled multiexponential signal to a 

5 transform operator constructed to provide siibstantially 

6 optimal linear resolution between outputs of the transform 

7 operator and an unknown model * 
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considered to go beyond the disclosure as filed (Rule 70.2(c)): 

(Any replacement sheet containing such amendments must be referred to under item 1 and annexed to this 
report.) 
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III. Non-establishment of opinion with regard to novelty, inventive step and industrial applicability 
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S the entire international application. 
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Re Item I 

Basis of the report 

1 . The basis of this international preliminary examination report is the application as 
originally filed. 

Re Item III 

Non-establishment of opinion with regard to novelty, inventive step and 
industrial applicability 

2. The present set of claims lacks clarity (Article 6 PCT) to an extent which does not 
permit an assessment of the invention of the basis of the subject-matter claimed. 
The reasons are as follows: 

2.1 The wording of the present independent claims can only be described as a 
juxtaposition of imprecise and vague expressions (.."a transform operator".. "an 
estimate of an unknown model". ."substantially optimal linear resolution"..) which 
does not convey any clear meaning to the skilled reader, because these 
expressions are used without any context and do not by themselves represent 
clear technical concepts. 

Claims should however be clear in themselves when being read with the normal 
skills including the knowledge about the prior art, but not including any knowledge 
derived from the description of the patent application - see PCT International 
Preliminary Examination Guidelines C-lll 4.2. 

2.2 Moreover, the claims also contravene the provisions of Article 6 PCT in that the 
presentation of 4 independent method claims (13,17,19 and 21) gives rise to two 
objections under Article 6 PCT, i.e. lack of conciseness and lack of clarity. 

2.3 As to conciseness, reference is made to the PCT Preliminary Examination 
Guidelines (PCT/GL Chapter III 4.1 .) in respect of the established practice that the 
requirement of conciseness applies not only to individual claims but to the claims 
as a whole. Rule 6.1 .a reinforces this conclusion. 
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2.4 The lack of clarity derives from the consideration that the prime function of the 
claims is to make clear what are the essential technical features of the matter for 
which protection is sought (cf. the first sentence of Art. 6 PCT). 
Present claims 13,17,19 and 21 appear in fact to provide four somewhat 
differently expressed versions of essentially the same (overly) broad features. 
These four alternative definitions leave the reader in doubt as to what are in fact 
the essential features and hence the primary purpose of Art. 6 PCT is not 
satisfied. 

3. In the present case, the plurality of independent claims drafted in vague and 
unspecific terms makes it impossible to establish an opinion with regard to 
novelty, inventive step or industrial applicability. 
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